home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 12 / CU Amiga Magazine's Super CD-ROM 12 (1997)(EMAP Images)(GB)[!][issue 1997-07].iso / CUCD / Sound / MusicIn / FFT_float_asm.a < prev    next >
Text File  |  1996-08-21  |  7KB  |  159 lines
  1. ;---------------------------------------------------------------------------
  2. ;
  3. ;       File:           FFT_float_asm.a
  4. ;
  5. ;       Author:         Stéphane TAVENARD
  6. ;
  7. ;       $VER:           FFT_float_asm.a  0.2  (20/08/1996)
  8. ;
  9. ;       (C) Copyright 1994-1996 Stéphane TAVENARD
  10. ;           All Rights Reserved
  11. ;
  12. ;       ----------------------------------------------------------------
  13. ;
  14. ;       #0 (23/04/1995) Initial revision
  15. ;       #1 (14/08/1995) Changed fmul by fsglmul & fdiv by fsgldiv
  16. ;       #1 (14/08/1995) Changed double COS,SIN by float COS,SIN
  17. ;       #2 (20/08/1996) Back to mul & div, cause '40 & '60 compat
  18. ;
  19. ;       Routines for FFT calculation with floating points values
  20. ;
  21. ;---------------------------------------------------------------------------
  22.  
  23.                   section  text,code
  24.  
  25.                   XDEF     @ASM_FFT_main_loop
  26.                   XDEF     _ASM_FFT_main_loop   ; for SAS-C 6.5
  27.                   XDEF     @ASM_FFT_energy_phi
  28.                   XDEF     _ASM_FFT_energy_phi  ; for SAS-C 6.5
  29. ;
  30. ;                 even
  31.  
  32. ;                 a0: float *x_real
  33. ;                 a1: float *x_imag
  34. ;                 a2: float *cos
  35. ;                 a3: float *sin
  36. ;                 d0: power
  37. ;
  38. _ASM_FFT_main_loop
  39. @ASM_FFT_main_loop
  40.                   movem.l  d2-d7/a2-a6,-(sp)
  41.                   move.l   a2,a4          ; a4 = cos
  42.                   move.l   a3,a5          ; a5 = sin
  43.                   moveq.l  #1,d1          ; d1 = d = 1
  44.                   clr.l    d2
  45.                   bset     d0,d2
  46.                   asr.l    #1,d2          ; d2 = li = 1<<( power-1 )
  47.  
  48.                   subq.w   #1,d0          ; d0 = i
  49. FFT_main_loop_1
  50.                   movem.l  a0-a1,-(sp)
  51.                   lea.l    (a0,d1.w*4),a2 ; a2 = &x_real[ d ]
  52.                   lea.l    (a1,d1.w*4),a3 ; a3 = &x_imag[ d ]
  53.  
  54.                   move.w   d2,d3          ;
  55.                   subq.w   #1,d3          ; d3 = j
  56. FFT_main_loop_2   move.w   d1,d4
  57.                   subq.w   #1,d4          ; d4 = l
  58. FFT_main_loop_3   fmove.s  (a0),fp0       ; fp0 = *real_ptr
  59.                   fmove.s  (a1),fp1       ; fp1 = *imag_ptr
  60.                   fmove.s  (a2),fp2       ; fp2 = *t_real_ptr
  61.                   fmove.s  (a3),fp3       ; fp3 = *t_imag_ptr
  62.                   fmove.x  fp2,fp4
  63.                   fmove.x  fp3,fp5
  64.                   fadd.x   fp0,fp2        ; fp2 = *real_ptr + *t_real_ptr
  65.                   fadd.x   fp1,fp3        ; fp3 = *imag_ptr + *t_imag_ptr
  66.                   fsub.x   fp4,fp0        ; fp0 = *real_ptr - *t_real_ptr
  67.                   fsub.x   fp5,fp1        ; fp1 = *imag_ptr - *t_imag_ptr
  68.                   fmove.s  fp2,(a0)+      ; fp2 -> *real_ptr++
  69.                   fmove.s  fp3,(a1)+      ; fp3 -> *imag_ptr++
  70.                   fmove.s  fp0,(a2)+      ; fp4 -> *t_real_ptr++
  71.                   fmove.s  fp1,(a3)+      ; fp5 -> *t_imag_ptr++
  72.                   dbra     d4,FFT_main_loop_3
  73.                   lea.l    (a0,d1.w*4),a0 ; real_ptr += d;
  74.                   lea.l    (a1,d1.w*4),a1 ; imag_ptr += d;
  75.                   lea.l    (a2,d1.w*4),a2 ; t_real_ptr += d;
  76.                   lea.l    (a3,d1.w*4),a3 ; t_imag_ptr += d;
  77.                   dbra     d3,FFT_main_loop_2
  78.  
  79.                   movem.l  (sp)+,a0-a1
  80.                   movem.l  a0-a1,-(sp)
  81.  
  82.                   asl.l    #1,d1          ; d = d<<1;
  83.                   asr.l    #1,d2          ; li = li>>1;
  84.                   beq      FFT_main_loop_9 ; li = 0 -> end of FFT
  85.                   lea.l    (a0,d1.w*4),a0 ; a0 = real_ptr = &x_real[ d ]
  86.                   lea.l    (a1,d1.w*4),a1 ; a1 = imag_ptr = &x_imag[ d ]
  87.                   move.w   d2,d3
  88.                   subq.w   #1,d3          ; d3 = j
  89. FFT_main_loop_5   move.l   a4,a2          ; a2 = cos_ptr = cos
  90.                   move.l   a5,a3          ; a3 = sin_ptr = sin
  91.                   move.w   d1,d4
  92.                   subq.w   #1,d4          ; d4 = l
  93. FFT_main_loop_6   fmove.s  (a0),fp0       ; fp0 = *real_ptr
  94.                   fmove.s  (a1),fp1       ; fp1 = *imag_ptr
  95.                   fmove.s  (a2),fp2       ; fp2 = *cos_ptr
  96.                   fmove.s  (a3),fp3       ; fp3 = *sin_ptr
  97.                   fmove.x  fp2,fp4
  98.                   fmove.x  fp3,fp5
  99.                   fmul.x   fp0,fp2        ; fp2 = (*real_ptr)*(*cos_ptr) (#2)
  100.                   fmul.x   fp1,fp3        ; fp3 = (*imag_ptr)*(*sin_ptr) (#2)
  101.                   fadd.x   fp2,fp3        ; fp3 = real
  102.                   fmul.x   fp1,fp4        ; fp4 = (*imag_ptr)*(*cos_ptr) (#2)
  103.                   fmul.x   fp0,fp5        ; fp5 = (*real_ptr)*(*sin_ptr) (#2)
  104.                   fsub.x   fp5,fp4        ; fp4 = imag
  105.                   fmove.s  fp3,(a0)+      ; *real_ptr++ = real
  106.                   fmove.s  fp4,(a1)+      ; *imag_ptr++ = imag
  107.                   lea.l    (a2,d2.w*4),a2 ; cos_ptr += li
  108.                   lea.l    (a3,d2.w*4),a3 ; sin_ptr += li
  109.                   dbra     d4,FFT_main_loop_6
  110.                   lea.l    (a0,d1.w*4),a0 ; real_ptr += d
  111.                   lea.l    (a1,d1.w*4),a1 ; imag_ptr += d
  112.                   dbra     d3,FFT_main_loop_5
  113. FFT_main_loop_9
  114.                   movem.l  (sp)+,a0-a1
  115.  
  116.                   dbra     d0,FFT_main_loop_1
  117.  
  118.                   movem.l  (sp)+,d2-d7/a2-a6
  119.                   rts
  120.  
  121. ;                 even
  122.  
  123. ;                 a0: float *x_real
  124. ;                 a1: float *x_imag
  125. ;                 a2: float *energy
  126. ;                 a3: float *phi
  127. ;                 d0: n
  128. ;
  129. _ASM_FFT_energy_phi
  130. @ASM_FFT_energy_phi
  131.                   subq.w   #1,d0
  132. FFT_energy_phi_0
  133.                   fmove.s  (a0)+,fp0
  134.                   fmove.s  (a1)+,fp1
  135.                   fmove.x  fp0,fp4
  136.                   fmove.x  fp1,fp5
  137.                   fmul.x   fp4,fp4           ; (#2)
  138.                   fmul.x   fp5,fp5           ; (#2)
  139.                   fadd.x   fp4,fp5           ; fp5 = en = real^2 + imag^2
  140.                   fmove.s  fp5,(a2)+         ; *energy++ = en
  141.                   ftst.x   fp5
  142.                   fbne.w   FFT_energy_phi_1  ; en <> 0
  143.                   fmove.s  #0,fp1            ; phi = 0
  144.                   bra.s    FFT_energy_phi_5
  145. FFT_energy_phi_1  ftst.x   fp0
  146.                   fbeq.w   FFT_energy_phi_3  ; real == 0
  147. FFT_energy_phi_2  fdiv.x   fp0,fp1           ; fp1 = imag / real (#2)
  148.                   fatan.x  fp1,fp1           ; fp1 = atan( imag / real )
  149.                   bra.s    FFT_energy_phi_5
  150. FFT_energy_phi_3  fmove.s  #1E-20,fp0
  151.                   bra.s    FFT_energy_phi_2
  152. FFT_energy_phi_5  fmove.s  fp1,(a3)+         ; *phi++ = phi
  153.                   dbra     d0,FFT_energy_phi_0
  154.                   rts
  155.  
  156. ;                 even
  157.  
  158.                   end
  159.